#ifndef INTERP_HPP
#define INTERP_HPP

#include <vector>
#include <stdexcept>
#include <iostream>

/**
 * @class NewtonInterp
 * @brief Implements Newton's interpolation formula.
 *
 * This class provides a way to compute interpolating polynomials for a given set of
 * data points (x, y) using Newton's divided differences approach.
 */
class NewtonInterp {
private:
    std::vector<double> x; ///< Vector of x values for the interpolation points
    std::vector<double> y; ///< Vector of corresponding y values
    std::vector<std::vector<double>> div; ///< Divided differences table used for computation

public:
    /**
     * @brief Constructor that initializes x and y values and computes the divided differences table.
     * 
     * @param xValues Vector of x values.
     * @param yValues Vector of y values corresponding to x.
     * 
     * Throws an exception if the sizes of xValues and yValues do not match.
     */
    NewtonInterp(const std::vector<double>& xValues, const std::vector<double>& yValues) {
        if (xValues.size() != yValues.size()) {
            throw std::invalid_argument("The size of x and y vectors must be the same.");
        }
        x = xValues;
        y = yValues;
        int n = x.size();
        div.resize(n, std::vector<double>(n, 0.0));
        computeDiv();
    }

    /**
     * @brief Compute the divided differences table based on x and y values.
     *
     * The divided differences table is filled using recursive relationships.
     * The first column is initialized with y values, and each subsequent column
     * is calculated based on the previous columns.
     */
    void computeDiv() {
        int n = x.size();
        // Initialize the first column with y values
        for (size_t i = 0; i < n; ++i) {
            div[i][0] = y[i];
        }
        // Calculate divided differences for each level
        for (size_t j = 1; j < n; ++j) {
            for (size_t i = 0; i < n - j; ++i) {
                div[i][j] = (div[i + 1][j - 1] - div[i][j - 1]) / (x[i + j] - x[i]);
            }
        }
    }

    /**
     * @brief Evaluate the Newton interpolation polynomial at a specified point.
     * 
     * @param point The x-coordinate at which to evaluate the interpolation polynomial.
     * @return The interpolated y-coordinate at the specified x.
     *
     * The evaluation uses the computed divided differences to generate the
     * interpolation polynomial value at the given point.
     */
    double evaluate(double point) const {
        int n = x.size();
        double result = div[0][0]; // Start with the first divided difference
        double product = 1.0;      // Holds the accumulated (point - x[k]) values
        for (size_t j = 0; j < n - 1; ++j) {
            product *= (point - x[j]);
            result += div[0][j + 1] * product;
        }
        return result;
    }
};

/**
 * @class HermiteInterp
 * @brief Implements Hermite's interpolation formula using both function values and derivatives.
 *
 * This class performs Hermite interpolation given values and derivatives at specified points.
 */
class HermiteInterp {
private:
    std::vector<double> x; ///< Vector of x values for interpolation points
    std::vector<double> y; ///< Vector of function values at x points
    std::vector<double> dy; ///< Vector of first derivatives at x points
    std::vector<std::vector<double>> div; ///< Divided differences table for Hermite interpolation

public:
    /**
     * @brief Constructor that initializes x, y, and dy values and computes the divided differences table.
     * 
     * @param xValues Vector of x values.
     * @param yValues Vector of y values corresponding to x.
     * @param dyValues Vector of first derivatives at x.
     * 
     * Throws an exception if the sizes of xValues, yValues, and dyValues do not match.
     */
    HermiteInterp(const std::vector<double>& xValues, const std::vector<double>& yValues, const std::vector<double>& dyValues) {
        if (xValues.size() != yValues.size() || xValues.size() != dyValues.size()) {
            throw std::invalid_argument("The size of x, y, and dy vectors must be the same.");
        }
        x = xValues;
        y = yValues;
        dy = dyValues;
        int n = 2 * x.size(); // Each x point has two entries in the table for Hermite interpolation
        div.resize(n, std::vector<double>(n, 0.0));
        computeDiv();
    }

    /**
     * @brief Compute the divided differences table for Hermite interpolation.
     *
     * Hermite interpolation requires the use of both values and derivatives at each point.
     * Each x point appears twice in the table, and special handling is done for the initial
     * divided differences using the derivative values.
     */
    void computeDiv() {
        int n = 2 * x.size();
        // Initialize the first column with function values and derivative information
        for (size_t i = 0; i < n; i+=2) {
            div[i][0] = y[i / 2];
            div[i + 1][0] = y[i / 2];
            div[i][1] = dy[i / 2];
        }
        for (size_t i = 0; i < n - 1; i+=2) {
            div[i + 1][1] = (y[i / 2 + 1] - y[i / 2]) / (x[i / 2 + 1] - x[i / 2]);
        }

        // Calculate higher-order divided differences
        for (size_t j = 2; j < n; ++j) {
            for (size_t i = 0; i < n - j; ++i) {
                div[i][j] = (div[i + 1][j - 1] - div[i][j - 1]) / (x[(i + j) / 2] - x[i / 2]);
            }
        }
    }

    /**
     * @brief Evaluate the Hermite interpolation polynomial at a specified point.
     * 
     * @param point The x-coordinate at which to evaluate the polynomial.
     * @return The interpolated y-coordinate at the specified x.
     *
     * This function uses the divided differences table to calculate the interpolation
     * polynomial value at the given point, considering both function values and derivatives.
     */
    double evaluate(double point) const {
        int n = 2 * x.size();
        double result = div[0][0]; // Start with the first divided difference
        double product = 1.0;      // Holds the accumulated (point - x[k]) values
        for (size_t j = 0; j < n - 1; ++j) {
            product *= (point - x[j / 2]);
            result += div[0][j + 1] * product;
        }
        return result;
    }
};

#endif // INTERP_HPP
